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ABSTRACT 

The morphology of the X-ray and radio emitting features in the central 50 kpc region 
around the galaxy M87 strongly suggests that buoyant bubbles of cosmic rays (inflated by an 
earlier nuclear active phase of the galaxy) rise through the cooling gas at roughly half the sound 
speed. In the absence of strong surface tension, initially spherical bubbles will transform into tori 
as they rise through an external medium. Such structures can be identified in the radio images of 
the halo of M87. During their rise, bubbles will uplift relatively cool X-ray emitting gas from the 
central regions of the cooling flow to larger distances. This gas is colder than the ambient gas and 
has a higher volume emissivity. As a result, rising "radio" bubbles may be trailed by elongated 
X-ray features as indeed is observed in M87. We performed simple hydrodynamic simulations to 
qualitatively illustrate the evolution of buoyant bubbles in the M87 environment. 

Subject headings: galaxies; active - galaxies; clusters; individual; Virgo - cooling flows - galaxies; 
individual; M87 - X-rays; galaxies 



1. Introduction 

The giant elliptical galaxy M87 continues to be 
the subject of numerous experimental and theo- 
retical studies. It is located at, or near, the center 
of the very nearby (~ 18 Mpc, 1 arcminute corre- 
sponds to ~5 kpc) irregular Virgo cluster and is 
at the center of a weak cooling flow with an es- 
timated mass deposition rate of some 10-40 Mq 
per year (e.g. Peres et al. 1998). The radio source 
Virgo A (3C274), associated with the galaxy M87, 
is well known because of the spectacular jet which 
is observed both in radio and optical bands (see 
e.g. Biretta 1999 for a review). In addition to the 
jet, which is contained within the central 2 kpc 
or 30" region, there is a lower surface brightness 
radio halo extending up to a distance of ~40-50 
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kpc from the nucleus (e.g. Kassim et al. 1993, 
Bohringer et al. 1995, Owen, Eilek and Kassim 
1999, 2000). 

The X-ray emission is strongly peaked at the 
position of the nucleus of M87 and is largely sym- 
metric around it. However, it was found in the 
Einstein HRI observations (Feigelson et al. 1987) 
and then in the ROSAT data (Bohringer at al. 
1995) that there are departures from a symmet- 
ric model and these departures resemble the mor- 
phology of some prominent features in the outer 
radio lobes. These findings initiated numerous dis- 
cussions on the interaction of the radio and X- 
ray emitting plasmas (e.g. Feigelson et al. 1987, 
Bohringer et al. 1995, Bohringer et al. 1999, 
Owen, Eilek and Kassim 1999,2000, Harris et al. 
1999). In a number of other cooling flow clusters 
with strong radio halos surrounding the dominant 
galaxy (e.g. Perseus - Bohringer at al. 1993, Mc- 
Namara et al. 1996, Churazov et al. 2000, Fabian 
et al. 2000a, A4059 - Huang and Sarazin 1998, 
Hydra A - McNamara et al. 2000) there is a clear 
anticorrelation of the radio and X-ray emitting 
plasmas. This makes the case of M87 especially 
interesting. 

We discuss below one particular model which is 
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able to reproduce qualitatively both the radio and 
X-ray morphologies observed in the central region 
surrounding the M87 galaxy. 

1.1. Radio and X ray morphology 

There is a long history of M87 radio halo ob- 
servations, starting with Mills (1952) shortly after 
the initial discovery of the radio source Virgo A 
(Bolton, Stanley and Slee 1949). Because of the 
very bright compact source, the dynamic range 
of the images was always a problem. Recently, 
Owen et al. (2000) presented a new high reso- 
lution, high dynamic range map of the halo of 
M87 at 327 MHz (sec Figure 1). For the sake 
of further discussion this map has been rotated 
by 90 degrees to have East up and North to the 
right. In this map several distinct structures are 
discernible. The high brightness centre represents 
the well-studied inner lobe structure (oriented ap- 
proximately east-west) with the famous jet point- 
ing west north-west (approximately bottom-right 
for the orientation of images adopted in Fig.l). 
The highly structured outer halo is much fainter. 
It consists of the ear-shaped (torus-like) eastern 
bubble, the much less well-defined western bubble, 
both of which are connected to the central emis- 
sion by a 'trunk', and the two very faint, almost 
circular, emission regions northeast and southwest 
of the centre. 

Owen et al. (2000) find a peak surface bright- 
ness in the region of the eastern bubble of about 
53 mjy/bcam while the circular regions range 
from 14 mJy/beam to about 25 mJy/beam. 
Rottmann et al. (1996) show that all of the 
described structures arc also present at 10.6 GHz. 
They also show that the outer halos have a radio 
spectral index in the range from —2.8 to —2.5 be- 
tween 4.8 GHz and 10.6 GHz. The spectral index 
is flatter by a factor of about 1.5 between 330 MHz 
and 1.5 GHz. This suggests that radiative energy 
losses of the synchrotron emitting particles lead 
to a significant steepening of the radio spectrum 
at about 3 GHz. Unfortunately, the maps used by 
Rottmann et al. (1996) have low resolution since 
they are obtained with single dish observations 
and so the spectral index distribution of the halo 
and its connection to the detailed structures, 
e.g., bubbles and attached "trunks", is not known. 

In Fig.l (lower left) we show the X-ray surface 



brightness distribution of the same region (and 
similar orientation) as the radio map. The data 
are the combined 200 ksec ROSAT/HRI image, 
adaptively smoothed to suppress high frequency 
noise. As pointed out by a number of authors 
(Feigelson et al. 1987, Bohringer at al. 1995, 
Owen, Eilek and Kassim 1999, Harris et al. 1999) 
there is evidence for a correlation between X-ray 
and radio emitting features. The simplest expla- 
nation of this correlation is that the excess X-ray 
emission is due to inverse Compton scattering of 
the cosmic microwave background photons by the 
same rclativistic electrons which produce the syn- 
chrotron radio emission (Feigelson et al. 1987). 
However, ROSAT/PSPC observations have shown 
that the excess emission has a thermal spectrum 
(Bohringer at al. 1995) and the X-ray emitting 
gas in these region has a lower temperature than 
that in the ambient regions. 

Thus, we summarize three observational facts 
which we are trying to explain below: 

• There are prominent "torus-like" features in 
the radio image. 

• There is a correlation (but not one to one) 
of the X-ray and radio bright regions. 

• Excess X-ray emission, associated with the 
radio bubbles and attached "trunks" , arises 
from thermal gas whose temperature is less 
than that of the ambient X-ray emitting gas. 

1.2. Analogy with atmospheric explosions 

Concentrating on the "torus-like" radio fea- 
tures one can note its striking similarity with some 
evolutionary stages of hot buoyant bubbles formed 
by a powerful (e.g. nuclear) explosion in the 
Earth's atmosphere. The transformation of the 
initially spherical bubble into a torus is a com- 
mon and well known property of buoyant bub- 
bles lacking strong surface tension (e.g. Walters 
and Davison 1963, Onufriev 1967, Zhidov et al. 
1977). Morphologically similar structures resem- 
bling " mushrooms" appear in Rayleigh- Taylor un- 
stable configurations: as the fluid rises through 
the ambient medium, Kelvin-Helmholtz instabili- 
ties create the torus-like head of the " mushroom" . 
Images of various laboratory experiments where 
such structures were observed and discussion of 
the relevant physical mechanisms are widespread 
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Fig. 1. Top-Left: 14'.6 x 16'. map of the radio halo of M87 at 327 MHz rotated 90°clockwise (North is 
to the right, East is up). The map was kindly provided by F.Owen (see Owen, Eilek and Kassim 2000 for 
original data). Top-Right: Adaptively smoothed ROSAT HRI X-ray image. The size and orientation of the 
image arc the same as for the radio image. Bottom: Possible geometry of the source inspired by analogy 
with "mushroom clouds" produced by powerful atmospheric explosions. The black region in the center 
denotes the inner radio lobes, gray "mushrooms" correspond to the buoyant bubbles already transformed 
into tori, and the gray lens-shaped structures are the 'pancakes' formed by the older bubbles. To explain 
the observed morphology, the source must be oriented close to the line of sight (dashed line). The pancakes 
are shown edge-on as shaded regions. 
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(see e.g. Batchelor 1967, Turner 1973, Inogainov 
1999). The idea that buoyancy plays an important 
role in the evolution of the radio lobes in galaxies 
was first proposed by Gull and Northovcr (1973) 
and has been used to estimate the life time of the 
radio lobes in M87 by Bohringer et al. (1995). 
The role of buoyancy was also discussed by Baum 
and O'Dea (1991) in application to PKS 0745-191. 
The similarity of the M87 "torus-like" features to 
the Raylcigh- Taylor mushrooms was pointed out 
by Churazov et al. (2000) and to a subsonic vortex 
ring by Owen, Eilek and Kassim (2000). 

Further pursuing the analogy with powerful ex- 
plosions, we note that during the transformation 
of a bubble into a torus some ambient gas is cap- 
tured and uplifted by the rising bubble/torus. In 
the cooling flow surrounding M87, the entropy of 
the X-ray emitting gas rises quickly with distance 
from the centre, and therefore the emissivity of 
the gas increases very steeply if it remains in 
pressure balance with the ambient gas (Bohringer 
et al. 1995, Nulsen 1997). This may qualitatively 
explain the correlation of the radio and X- 
ray emitting plasmas and naturally accounts for 
the thermal nature of the excess thermal emission. 
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Fig. 2. — Initial distribution of electron den- 
sity, temperature, pressure and gravitating mass 
(adopted from Nulsen and Bohringer 1995) as- 
sumed in the simulations. 



In the case of an atmospheric explosion, the fi- 
nal stage of the hot bubble occurs when the bubble 
reaches a height at which the density of the am- 
bient gas is equal to the density of the bubble. 
Then (neglecting oscillations near the equilibrium 
position) the bubble expands laterally (forming a 
"pancake") to occupy a thin layer in the atmo- 
sphere having the same density. The bubbles/tori 
in M87 may share the same fate. In a spheri- 
cally symmetric potential, the bubble will try to 
fill a segment (equipotential surface) of a sphere. 
The largest distinct features in the radio map (see 
Fig.l) could be just those late-stage bubbles. 

The sketch of the possible overall source struc- 
ture based on the analogy with a powerful atmo- 
spheric explosion is shown in Fig.l. Here the black 
region in the center denotes the inner radio lobes, 
the gray mushrooms correspond to the buoyant 
bubbles already transformed into tori and the gray 
lens-shaped structures are the "pancakes" formed 
by the older bubbles. 

Buoyant bubbles and their role in the complex 
X-ray and radio morphology of M87 have been dis- 
cussed previously e.g. in Bohringer et al. (1995), 



Nulsen (1997), Bohringer (1999), Owen, Eilek and 
Kassim (2000). Below wc follow the line of argu- 
ments suggested in Churazov et al. (2000) and 
present hydrodynamic calculations of the radio 
lobe evolution to illustrate the qualitative picture 
described above. 

2. Method 

In the next sections we describe the method, 
initial conditions, and assumptions used in the nu- 
merical simulations of bubbles in a hot gaseous at- 
mosphere. We first detail the method and the use 
of "tracer" particles to track the motion of hot gas 
which is entrained in the rising bubble. We next 
describe our approach to calculation of the syn- 
chrotron radio emission from the evolving bubble. 

2.1. Hydrodynamic simulations 

The simulations were obtained using the 
ZEUS-3D code which was developed especially 
for problems in astrophysical hydrodynamics 
(Stone & Norman 1992a). The code uses finite 
diff'erencing on an Eulerian grid and is fully 
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explicit in time. It is based on an operator-split 
scheme with piecewise linear functions for the 
fundamental variables. The fluid is advected 
through a mesh using the upwind, monotonic 
interpolation scheme of van Leer (1977). For a 
detailed description of the algorithms and their 
numerical implementation sec Stone & Norman 
(1992a, b). To study the spectral aging of the 
relativistic gas, the ZEUS code was modified to 
follow the motion of 'tracer' particles which are 
advected with the fluid. 

In our simulations we employed an ideal gas 
equation of state and we ignored the effects of 
magnetic flelds and rotation. The cooling time 
is sufHciently long compared to the time scales 
considered here since, for an electron density of 
0.02 cm~^ and temperature of the order of 2 keV 
(see Fig. 2) the cooling time is ~ 5 x 10*^ years. 
This is an order of magnitude longer than the 
typical duration of the run of ~ 5 x 10*^ years. 
We therefore completely neglected cooling. The 
simulations were computed on a spherical grid in 
two dimensions and were performed on an IBM 
RS/6000 cluster. The computational domain 
spans 51 kpc in radius and 7r/2 rad in angle and 
was covered by 200 x 200 grid points (in the r- 
and ^-direction). 

A model for the mass and initial temperature 
was adopted from Nulsen & Bohringer (1995) 
which is shown in Fig. 2. The gas density distri- 
bution was then found by assuming hydrostatic 
equilibrium to maintain an initially static model. 
A spherical bubble was created with a radius of 
Tb = 5 kpc at a distance of d = 9 kpc from the 
gravitational center. It was made buoyant by re- 
ducing its density with respect to the background 
by a factor of 100 and simultaneously raising 
the temperature by the same factor. During 
its subsequent evolution, the gas is treated as a 
single fluid which is assumed to obey a polytropic 
equation of state with F = 5/3. The bubble 
was flUed with 2000 tracer particles and was 
placed on the 9 = axis along which it rose 
in the gravitational potential. Subsequently, we 
rotated the two-dimensional output about this 
axis (azimuthally). To test the assumption of 
rotational symmetry, we performed some 3D 
simulations with reduced resolution and found 



that this assumption is fairly robust as far as the 
global motions and morphologies are concerned. 
As we discuss in detail below, buoyancy deforms 
the bubble and drives it through the ambient 
medium as shown in Fig. 3. Shear instabilities 
cause the formation of tori which separate from 
the main bubble. This process, also called 
"vortex shedding" , has been observed and studied 
extensively (see Norman et al. 1982). 

Finally, we should address some issues related 
to the accuracy of these kinds of finite-difference 
hydrodynamic simulations. While the code can 
simulate large-scale mixing due to Rayleigh- Taylor 
and Kelvin-Helmholtz instabilities, it does not in- 
clude real particle diffusion. Any observed diffu- 
sion is therefore entirely numerical. The boundary 
between the bubble and the ambient medium be- 
comes less sharp as the simulation proceeds due 
to discretization errors in the advection scheme. 
For a test of the advection algorithm in the ZEUS 
code see Stone & Norman (1992a, b). In simple 
advection tests, it was found that, during the ad- 
vection of a sharp discontinuity over a grid of 200 
zones, the discontinuity is spread over 3-4 grid 
cells. Therefore, the small features in our simu- 
lation are likely to be affected by these advection 
errors whereas the larger features are not. 

Second, numerical viscosity is also responsible 
for suppressing small-scale instabilities at the in- 
terface between the bubble and the cooler, sur- 
rounding. X-ray emitting gas. To assess the effects 
of numerical viscosity, we have repeated our simu- 
lations on grids with 100 x 100 grid points. From 
our experiments, we can conclude that "global pa- 
rameters" such as the position and size of the bub- 
ble as well as the presence of "toroidal" structure 
are relatively insensitive to the resolution. The de- 
tailed small scale morphology does depend on the 
resolution and the initial conditions. 

Finally, we stress that the simulations presented 
here are not intended to provide a detailed sim- 
ulation of the radio/X-ray structure surrounding 
M87. Rather they are intended as a guide, a toy 
model, to assist our interpretation of the complex 
emission (see discussion in Section 3.1). 
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2.2. Simulation of the radio emission 

The observations of radio synchrotron emission 
from the buoyant bubbles in M87 imply the pres- 
ence of a relativistic, magnetized plasma. How- 
ever, the simulations presented in previous sec- 
tions are based on a purely hydrodynamic scheme 
with a single, non-relativistic fluid. To estimate 
the radio flux expected from the simulation, we 
make the following assumptions. 

The relativistic plasma and the associated mag- 
netic field, the 'relativistic fluid', are confined 
within small volumina intermixed with the non- 
relativistic, thermal plasma governing the dynam- 
ics. This implies that the relativistic fluid may 
mix with the thermal fluid on macroscopic scales 
but not on microscopic ones. This is very simi- 
lar to the behavior of air bubbles in water which 
are confined by surface tension (e.g. Zhidov et al. 
1977). 

We assume the magnetic field to be tangled on 
scales small compared to the size of the volumina 
of relativistic fluid. This allows us to treat the 
magnetic field as part of the relativistic fluid. The 
total energy density inside the small volumina is 
therefore the sum of the energy density of the mag- 
netic field, Mb, and that of the relativistic parti- 
cles, Uc- We assume that the small volumina of 
relativistic fluid are in pressure equilibrium with 
the thermal fluid throughout the simulation, i.e. 
3pth = ub + Uc- We also assume that there is no 
significant re-acceleration of relativistic particles 
in the bubble during its evolution. 

At the start of the simulation, we assume some 
initial value of the magnetic field (the same value 
for every small volumina of relativistic particles 
within the bubble) and a power law distribution 
for the relativistic electrons (in terms of the par- 
ticle Lorentz factor), 7, ne{to)d'y = no'y~^d'y. At 
the start of the simulation is therefore given by 



mf.(?n, 



/■Tmax 
•'Tmin 



7-^(7-l)d7, (1) 



where wc assume 7,nin = 1, 7max = 10®. For a 
given initial magnetic field the normalization Uq 
of the particle distribution is then adjusted such 
that the entire relativistic fluid is in pressure equi- 
librium with the thermal fluid. 

For values oip > 2, the energy density of rela- 
tivistic particles is dominated by low energy par- 



ticles. The life time of such particles, due to ra- 
diative losses, is much longer than the duration of 
the simulation (and no Coulomb losses are present 
since we assume that thermal gas is not mixed 
with relativistic particles on microscopic scales). 
The motion of the bubble (see Section 3) is found 
to be subsonic even with respect to the cluster 
gas. We therefore assume that the energy den- 
sity of the magnetic filed and relativistic parti- 
cles changes adiabatically according to the changes 
of the pressure of the surrounding thermal gas: 
Mb oc Me oc pth- Note, however, that our hydro- 
dynamic simulations assume an adiabatic index 
of 5/3 for the nonrelativistic gas, while the rel- 
ativistic fluid expands and contracts adiabatically 
with an adiabatic index = 4/3. This implies 
that, during the expansion of a buoyant bubble, 
the fraction of the volume occupied by the rela- 
tivistic fluid increases relative to the volume of the 
thermal fluid as 14 oc V^^^. Strictly speaking this 
means that accounting for relativistic fluid with a 
different adiabatic index will not affect the dynam- 
ics of the bubble evolution as long as the fraction 
of the volume occupied by the relativistic fluid is 
small. We will sec in the following that this as- 
sumption is probably invalid, but we believe that 
it does not very severely affect the results. 

The energy losses of the relativistic particles 
are due to adiabatic expansion, synchrotron radia- 
tion and inverse Compton scattering of the CMB. 
These losses depend on the pressure history of each 
volumina of relativistic particles. The pressure 
history can be recovered using the tracer particles. 
For a given region within the bubble the change of 
the Lorentz factor 7 of a given electron due to adi- 
abatic expansion or contraction of the relativistic 
fluid and radiative losses can be written as 



^7 



7 dpt\,{t) 



Pth(i) 



3meC 



7 



dt 



WB(io) 



Pth(^) 



■ UlC 



(2) 



wherea=(r-l)/r = l/4. 

This can be integrated (e.g. Kaiser, Dennett- 
Thorpe & Alexander 1997) to give the relation be- 
tween the Lorentz factor of the electron at time t 
and the initial Lorentz factor of the same electron 
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with 



b = 



4crT 

3mpC 



MB(io) 



Pth(to) 



UlC 



pth(t)"rfi. (4) 



The vahio of b is calculated from the pressure his- 
tory of the tracer particles. The integral in this ex- 
pression must be evaluated numerically using the 
discrete time steps of the simulation. 

Finally, the adiabatic change of the volume in 
which the relativistic particles are contained im- 
plies Uoit) = no{to)\pth{t)/pth{to)V^^ ■ The en- 
ergy distribution of relativistic electrons at time t 
is then 



Pth{to)T 

(5) 

The energy distribution of the relativistic par- 
ticles combined with the strength of the magnetic 
field yields the monochromatic synchrotron emis- 
sivity (e.g. Longair 1994) 



Jv = -^cttcub 



7^(/)(7, 1/) ne d7, (6) 



where ^(7, p) is the dimensionless spectrum emit- 
ted by a single electron. 

To obtain the synchrotron surface brightness 
from the simulation, we assume rotational sym- 
metry about the polar-axis and integrate the emis- 
sivity along lines of sight through the resulting 3- 
dimensional model. 

3. Results and discussion 

With the foundation and assumptions de- 
scribed in the previous section, we describe the 
results of two numerical simulations. Wc; first 
present the results from the two numerical simula- 
tions with different initial conditions for the bub- 
bles. We compute the velocity and position versus 
time of the rising bubble. For the first simulation, 
we compute the synchrotron radio emission from 
the bubble for two orientations, one aligned with 



the symmetry axis in the plane of the sky and the 
second withthe symmetry axis aligned 45° with re- 
spect to the plane of the sky. The two orientations 
differ in the initial conditions for the magnetic field 
strength in order to reproduce the observed ra- 
dio spectrum. We derive the X-ray appearance of 
the bubble from the first simulation, also for the 
two different orientations. Finally, we discuss the 
"pancake" phase and calculate the largest distance 
to which bubbles can rise as a function of their ini- 
tial location and volume fraction of trapped, hot 
gas. 

3.1. Bubble configuration 

Fig. 3 shows snapshots of the temperature 
distribution during the evolution of the bubble. 
One can identify two stages of the bubble evolu- 
tion. The initially spherical bubble fiattens and 
develops a "cone" at the rear, which is filled with 
entrained gas. The bubble then transforms into 
a torus which at later stages may fragment into 
smaller structures, but which continue to show a 
global toroidal structure. Note that ambient gas 
(in particular the gas captured during the bubble 
— > torus transformation) occupies the central part 
of the rising structure. Note also that in the last 
stages (shown in Fig. 3) the coldest temperatures 
are found not at the center of the cluster but 
within the rising structure. We stress again 
that no radiative cooling was included in the 
simulations and these low temperatures regions 
are associated with gas which has been uplifted 
from the central regions and has expanded adi- 
abatically to match the ambient gas pressure at 
the current location of the bubble. 

Of course the particular shape of the 
torus/mushroom structure (e.g. two tori at dif- 
ferent distances from the center in Fig. 3) is the 
result of our choice of initial and boundary condi- 
tions and assumed axial symmetry. For example, 
Fig. 4 shows snapshots of the temperature distribu- 
tion for slightly different initial conditions (with an 
initial temperature contrast, between the ambient 
gas and the bubble, which is a factor of two larger 
than in the previous run). Note the changes in the 
detailed morphology of the rising bubble compared 
to Fig. 3. In general one would expect a large va- 
riety of structures to be formed, especially in 3D. 
But the torus-like geometry is a generic feature 



7 





/ \ 



Fig. 3. — Temperature distribution in the gas at 5 time steps. From left to right: 0, 8.4, 21, 42 and 67 Myrs 
after the start of the simulation. Each box is 40 by 20 kpc. The center of the cluster is at the bottom of the 
box. Temperature is color coded in the 0.7 (black) to 5 keV (white) range. The temperature of the cluster 
thermal gas (blue) changes from ^ 1 keV at the center to ~ 1.7 keV at a distance of 40 kpc. All temperatures 
above 5 keV are white. Thus, the hot "radio-emitting plasma" which initially has a temperature of order 100 
keV is white. Numerical diffusion (due to very strong gradients in density and temperature) creates regions 
with intermediate values of temperature (yellow, red, green) along the boundaries of the bubble. Note that 
during later stages, the coldest gas is not at the centre of the cooling flow, but is associated with the rising 
bubble. Since radiative cooling was not included in the simulations these cold regions represent uplifted (and 
adiabatically expanded) gas. 



of rising bubbles in the absence of strong surface 
tension. 

3.2. The velocity of the rising bubble 

The simplest velocity estimate of the rising bub- 
ble can be obtained by equating the ram pressure 
and buoyancy forces acting upon a bubble (e.g. 
Gull and Northover 1973). The buoyancy force is 
obviously: 

Fy, = Vg{p,-pi,), (7) 

where V is the bubble volume, g is the gravita- 
tional acceleration (we assume the ambient gas is 
in hydrostatic equilibrium), pa and pb are the mass 
densities of the ambient and the bubble gas respec- 
tively. The ram pressure (drag) is: 

Fd ~ C^Sv^p,, (8) 

where S is the cross section of the bubble. The 
numerical coefficient C (drag coefficient) depends 
on the geometry of the bubble and the Reynolds 



number. Thus the terminal velocity of the bubble 
is 

/ 2 Pa - Pb 

Here the factor of (pa — Pb)/pa can be dropped if 
the bubble density is low compared to the ambi- 
ent gas density. The expression for the terminal 
velocity can be further rewritten using the Kep- 
lerian velocity at a given distance from the clus- 
ter center: v ~ ^ {r/R){8/3C)vK, where r is the 
bubble radius, R is the distance from the center 
and vk = VgR is the Keplerian velocity. For 
a solid sphere moving through an incompressible 
fluid the drag coefficient C is of the order 0.4-0.5, 
for Reynolds numbers in the range ~ 10^-10^ (e.g. 
Landau and Lifshitz 1963). The drag coefficient 
for a rising bubble should of course be different. 
First of all, it is not a solid sphere, but a rather 
complicated structure, resembling a smoke ring in 
air. Such rings may travel large distances through 
the air with low drag. Secondly the Mach number 
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Fig. 4. — The same as in the previous figure but with an initial temperature contrast, between the ambient 
gas and the bubble, which is a factor of two larger. Note that detailed morphology of the rising bubble is 
rather sensitive to the initial conditions (and boundary conditions, symmetry assumed etc. ) and can be 
very different from what is shown in these figures. The presence of the torus-like structure is, however, a 
generic feature of the rising bubble phenomenon. 



of the rising bubble in our simulations is ^^0.6-0.7 
and the compressibility of the gas is important. 
This would tend to increase the drag coefficient, 
e.g., for a sphere moving at M ^ 0.7 the drag coef- 
ficient is C ~ 0.6. Thirdly, in a stratified medium 
an additional contribution to the drag may come 
from the excitation of internal gravity waves. In 
our simulations, the typical Keplerian velocity was 
~ 400 km s^^. Fig. 5 shows the position of the 
bubble front as a function of time. The solid line 
shows the expected position of the feature moving 
with the constant velocity of ~ 390 km s~^. The 
effective drag coefficient which can be estimated 
from the simulations, assuming the above formula 
for the rise velocity, is C ^ 0.75. We note here 
that this value may in turn be affected by the nu- 
merical resolution adopted here. However, since 
the velocity of the bubble depends on the square 
root of the drag coefficient, formula (9) still yields 
crude order of magnitude estimates of the bubble 
velocity. Thus, a large and strongly underdense 
bubble will rise with a velocity comparable to the 
Keplerian velocity. 

In the case of M87, the bubble's environment 
may not be very underdense compared to the am- 
bient gas (since no clear X-ray holes have been 
seen) and the velocity of the bubble may be some- 
what smaller. However, we believe that the above 



estimate of the velocity is approximately correct 
since the velocity depends only weakly (as the 0.5 
power) on the bubble size and gas density gradi- 
ents. This velocity defines a time scale of ~ few 
10^ years for the evolution of the bubble to the 
stage seen as a bright torus in Fig.l. 

3.3. The rising bubble in synchrotron 
emission 

In Figure 6 we show the simulation results 
at 327 MHz. The radio emission was calculated 
for the hydrodynamic simulations shown in 
Fig. 3. Upper and lower rows correspond to 
different viewing angles. The symmetry axis is 
perpendicular to the line of sight in the upper 
row (case 1) in this figure and at an angle of 45° 
to the line of sight (case 2) in the lower row. 
From the radio observations of M87 we note that 
the well-defined eastern bubble is observed at a 
distance of roughly 25 kpc from the source centre. 
This corresponds to the third panel in the upper 
row of Figure 6, i.e. bubble age of 42 Myrs, while 
for the simulation projected at an angle of 45° 
(lower row), the fourth panel, at an age of 67 
Myrs, gives the correct projected distance. Note 
here that these distances are for the primary 
bubble/ vortex ring. Figure 6 clearly shows the 
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Fig. 5. — Position of the bubble front as a func- 
tion of time. For comparison solid line shows the 
motion with the constant velocity of 390 km/s. 

presence of a secondary vortex ring which is closer 
to the centre. The appearance of this secondary 
ring is fairly dependent on the particular choice of 
initial conditions (see Section 2 and Figures 3,4). 
In the case studied here, the synchrotron surface 
brightness of the secondary vortex ring exceeds 
that of the primary bubble. No such structure 
is observed in M87 and so we concentrate in the 
following discussion on the primary bubble. 

For an initial power law distribution of relativis- 
tic electrons in a constant magnetic field, a signif- 
icant steepening of the synchrotron spectrum will 
be observed for frequencies higher than the break 
frequency, inbreak, which is given by (e.g. Leahy 
1991) 

i/break = 2.52 X 10^— GHz, (10) 

(S' + ^cmb) 

where B is the strength of the magnetic field in 
/uG and t is the time (in Myrs) for which the 
plasma was exposed to this field. Bqmb is the 
strength of the magnetic field in /zG equivalent to 
the energy density of the CMB, uic- For the red- 
shift of M87, uic ~ 4.2 X 10"^^ erg cm'^ and so 



-BcMB ^ 3.2/xG. From equation 10 the time t can 
be expressed as a function of the break frequency 
and the strength of the magnetic field. For any 
given break frequency, the time t is maximal if 
-B(tmax) = Bqmb/V^ ~ 2;uG. We have seen that 
the observed radio spectral indices (Section 1.1) 
require inbreak ^ 3 GHz, for which a maximum age 
is tmax ~ 10* years. This is comparable with the 
time span of our simulations. For example, if the 
symmetry axis is at angle of 45° with the line of 
sight (case 2 above), then the front of the bubble 
will reach a projected distance of ~25 kpc after 
67 Myrs. This simple estimate shows that we re- 
quire a rather low value of the initial magnetic 
field strength (comparable to B(tmax) ^ 2/iG) in 
order to ^void having the i^reak 

lie at a frequency 

lower than ~ 3 GHz. 

For case 1, we used p = 2.3 and an initial 
magnetic field strength of 6.5 fiG. We then cal- 
culated the synchrotron emission of the bubble as 
described in Section 2.2. After ^ 42 Myrs, i.e. 
when the bubble is ^25 kpc away from the cen- 
ter, the spectral index of the primary bubble is 
-2.7 from 4.7 - 10.6 GHz and -1.0 from 330 - 
1500 MHz which compares well with the observed 
values (Rottmann et al. 1996, Section 1.1). 

For case 2 the magnetic field inside the bubble 
must be lower since here the bubble is older when 
its upper edge reaches the correct projected height 
from the cluster centre. We set B = 2.9/iG and 
p = 2.3. For this case the high frequency spectral 
index is found to be —3.0 while at lower frequen- 
cies it is —1.0. 

The assumption of pressure equilibrium deter- 
mines the initial energy density of the relativistic 
particles. The energy density of the magnetic field 
is 1.7 X 10~^^ ergs cm~^ for case 1 and 3.4 x 10~^^ 
ergs cm~^ for case 2. This is much lower than 
the energy density of the thermal cluster gas 
and the particles energy density is therefore 
tic ^ 3pth ^ 4.2 X 10"^^° ergs cm"'' for both 
cases discussed. This implies a severe departure 
from the equipartition condition of Ub ~ We often 
assumed for synchrotron emitting plasmas. The 
initial volume filling factor of the radio emitting 
plasma in the initial bubble can then be adjusted 
to provide the required surface brightness of 
the final bubble. The peak surface brightness 
predicted in the region of the primary vortex ring 
is 55 mJy/pixel in case 1 (see Owen et al. 2000, 
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Fig. 6. — Synchrotron surface brightness at 327 MHz at four time steps for the simulation shown in Fig. 3. The 
upper row (case 1) shows the buoyant bubble with a symmetry axis in the plane of the sky and B ~ 6.5/iG 
initially. The lower row (case 2) is for a symmetry axis inclined 45° to the line of sight and B ~ 2.9y^G 
initially. From left to right both rows show the bubble at 4.2, 21, 42 and 67 Myrs after the start of the 
simulation. The filled contours bound regions of 0.017 (black), 0.04, 0.09, 0.2, 0.45, 1.0, 2.4, 5.5, 13, 28 and 
66 (white) mJy per pixel, where a pixel is 250 x 250 pc at the distance of M87 (~18 Mpc). 



Section 1.1), however, the required volume filling 
factor of the relativistic plasma is very close to 
unity. For case 2, even for a filling factor of unity, 
the peak surface brightness of the bubble is only 
2.4 mjy/pixel. Furthermore, we neglected any 
possible contribution from the relativistic ions 
(the so-called "k" factor) and thermal plasma 
to the pressure inside the bubble, which would 
further reduce its synchrotron emissivity. We 
also note here that volume filling factors ^1 pose 
a serious problem for our assumption that the 
relativistic fluid does not influence the overall 
dynamics of the buoyant bubble. 

The problems discussed above arise from the life 
time of the synchrotron emitting particles being 
short compared to the dynamical time scale of the 
structure. Gull & Northover (1973) first noted 
this lifetime problem for buoyant models. Below, 
we discuss various solutions. 

From Figure 6, we find that the diameter of 



our primary bubble is smaller than the observed 
diameter of the eastern bubble in M87. As we 
showed in Section 3.2, the rise time of the bub- 
ble is correlated with its size. It is therefore quite 
possible that the eastern bubble in M87 is some- 
what younger than our simulations suggest. Such 
a younger age would allow a stronger magnetic 
field in the bubble and could explain the observed 
spectrum, which in turn would allow a smaller vol- 
ume filling factor for the relativistic fluid. 

The buoyant bubble may also be filled mainly 
by a weak magnetic field interspersed with regions 
of a much higher field strength but small volume 
filling factor. The relativistic electrons may then 
survive for a much longer time in the mainly weak 
magnetic field and only diffuse slowly into the re- 
gions of greater field strengths which are responsi- 
ble for most of the observed synchrotron emission 
(Eilek, Melrose & Walker 1997). Subsonic turbu- 
lence can also enhance the energy density of the 
magnetic field in regions with initially weak field 
(e.g. Eilek, Owen and Zhou 1999). 
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If the magnetic field in the buoyant bubbles is 
stronger than that estimated above and rather ho- 
mogeneous, then in situ acceleration of relativis- 
tic particles is required. In this case, the age of 
the buoyant bubble cannot be constrained from ra- 
dio observations. The two latter mechanisms were 
also discussed by Owen ct al. (2000) to explain 
the radio emission of the outer M87 halo. 

The fundamental conclusion of these simula- 
tions is that one can reproduce the observed ra- 
dio brightness and the break in the spectrum of 
the torus-like feature in M87 assuming that (i) 
the bubble, filled with the relativistic plasma, is 
in pressure equilibrium with the thermal gas, (ii) 
the initial magnetic field strength is relatively low 
so as to prolong the lifetime of the electrons and 
(iii) no in situ reacceleration of particles or gen- 
eration of magnetic field occurs during the bubble 
evolution. However, the equipartition condition is 
sacrificed and the parameter space which satisfies 
these assumptions is rather tight. It is likely that 
the real situation is much more complex than was 
assumed in the simple model discussed above. 

3.4. The rising bubble in X— ray emission 

As is clear from Fig. 3, the bubble captures am- 
bient gas during its transformation to a torus and 
then uplifts this gas to larger radii. Since all mo- 
tions in our simulations are subsonic, this uplifted 
gas remains in approximate pressure equilibrium 
with the ambient gas at larger radii. The volume 
emissivity of this gas can be compared with the 
volume emissivity of the ambient gas: 



cu = n^A(T„) 
ca = nlA{Ta), 



(11) 



where n and T arc the density and temperature 
of the uplifted and ambient gas respectively, A(T) 
is the emissivity in a given energy band. If we 
assume adiabatic evolution of the entrained and 
uplifted gas then: 
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Fig. 7. — Volume emissivity of the uplifted and 
adiabatically expanded gas relative to the ambi- 
ent gas. The curve is normalized to unity at 1 
kpc. The enhancement factor, F(r), is the ratio 
of the emissivity of the gas uplifted from 1 kpc 
to r to that of the ambient gas at r. The ra- 
tio F{r2)/ F{r\) characterizes the factor by which 
volume emissivity of the gas uplifted from ri to r2 
exceeds the volume emissivity of the ambient gas 
at r2- 

where the subscript stands for the values of the 
uplifted gas at its initial location where it was cap- 
tured by the rising bubble. Thus the uplifted gas 
will have an emissivity which is higher than the 
ambient gas by a factor: 
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The enhancement factor F (excluding the cooling 
functions) is shown in Fig. 7. One can see that 
in a cooling flow (density decreases with distance 
while temperature increases), the emissivity of the 
uplifted gas greatly exceeds the emissivity of the 
ambient gas if the displacement is large. The ratio 
A(r„)/A(ra) (when the temperatures are near 1 
keV and the soft X-ray band is considered) further 
increases the brightness of the uplifted gas with 
respect to the ambient gas. 
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Fig. 8. — Expected X-ray morphology of the bubble for the same times as in Fig. 6. Each box is 40 by 40 kpc. 
The figure shows the relative deviation of the X-ray surface brightness (in the ROSAT energy band) with 
respect to the unperturbed X-ray emission of the cooling flow. The darker regions are X-ray underluminous 
and lighter regions are X-ray bright. Black corresponds to regions dimmer by 40% (or more) than the 
unperturbed value and white corresponds to regions brighter by 40% (or more) than the unperturbed value. 



To illustrate how the X-ray surface brightness 
is affected by the a radio emitting plasma bub- 
ble and lumps of uplifted cold gas, we calculated 
the projected surface brightness distribution (for 
the ROSAT energy band from 0.5 to 2 keV) for 
several stages of the bubble evolution. When cal- 
culating these maps, we assumed axial symmetry 
and assumed that the axis of the bubble is either 
perpendicular to the line of sight (upper row of im- 
ages in Fig. 8) or that the angle between the line 
of sight and bubble axis is 45° (lower row of im- 
ages in Fig. 8). Fig. 8 shows the relative deviation 
of the X-ray surface brightness from the unper- 
turbed value of the surface brightness at the same 
distance from the center. As expected the "trunk" 
of the mushroom is brighter than the surrounding 
regions. 

Fig. 9 illustrates the effect of entrainment of 
cold gas as the bubble rises. Shown in this fig- 
ure arc the distributions of tracer particles at the 
initial and final (after ^67 Myrs) stages of the 
bubble evolution. The tracer particles fill the en- 
tire volume outside the bubble. The color assigned 



to each particle characterizes its initial distance 
from the cluster center. For example, "green" and 
"blue" particles seen in the upper torus have been 
uplifted from the central region of the galaxy or 
cluster atmosphere. 

Abundance gradients are frequently observed in 
cluster atmospheres around central bright galaxies 
(e.g. Matsumoto et al. 1996). If our interpretation 
of the X-ray features as uplifted lumps of gas is 
correct, then one would expect the abundance in 
these regions to be higher than in the surrounding 
regions. This can be tested with Chandra and 
XMM. 

3.5. The pancake stage 

If the lumps of radio and X-ray emitting plas- 
mas in the bubble are bound together (e.g. by 
magnetic stresses as in the multiphase cooling flow 
picture of Nulsen (1986)) then one can calculate 
the maximum radius to which the whole bubble 
can rise for given volume fractions of "entrained" 
gas and (weightless) radio emitting plasma (see 
eq. (6) in Churazov et al. 2000). Fig. 10 shows 
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Fig. 9. — Initial (top) and final (bottom) positions of tracer particles, originating within different regions 
of the cluster gas. Each particle is color coded according to its initial distance from the center. Bottom 
panel shows the positions of particles ~67 Myrs after the beginning of the simulation. "Green" and "blue" 
particles seen in the upper torus (bottom panel) have been uplifted from the central region of the cooling 
flow. 



the dependence of the final position of the bub- 
ble on the initial position and fraction of the bub- 
ble volume occupied by entrained gas which was 
captured near the initial location of the bubble. 
At this final location, the average mass density of 
the bubble (consisting of separate lumps of radio 
and X-ray emitting plasma) is equal to the den- 
sity of the ambient gas. At the final position the 
bubble will spread along the equipotential (isoden- 
sity) surface and may form a pancake-like struc- 
ture. The X-ray appearance of the bubble at this 
stage crucially depends on whether mixing of ra- 
dio and X-ray emitting plasmas is microscopic or 
macroscopic (Bohringer et al. 1995, Nulsen 1997) 
A simple simulation demonstrating pancake 
formation is presented in Fig.ll. In this simula- 
tion, we assumed that the bubble is initially filled 



with gas having the same entropy as the cluster 
gas at a distance of ~25 kpc. This is done by set- 
ting the temperature of the gas inside the bubble 
a factor of ~2 higher than the temperature of the 
surrounding cluster gas. The density in the bub- 
ble was set a factor of 2 lower in order to main- 
tain pressure equilibrium. The resulting modest 
contrast in entropy between the bubble and the 
surrounding gas approximately models the situa- 
tion when the filling factors of the radio emitting 
plasma and thermal plasma in the bubble are com- 
parable as considered in the previous section. The 
upper row of Fig.ll shows the entropy distribu- 
tion at the beginning of the simulation and after 
~ 2.6 X 10^ years. In the lower panel, we show the 
distribution of tracer particles, which initially are 
located within the bubble, at the same times. The 



14 




Initial Radius, kpc 

Fig. 10. — The maximum radius to which the bub- 
ble can rise depending on its initial location and 
the volume fraction of gas, captured during the; 
transformation of the bubble into a torus. It is as- 
sumed that the bubble consists of separate lumps 
of radio and X ray emitting plasmas bound to- 
gether. Each curve corresponds to a given volume 
fraction (labels under the curves) of entrained gas 
added to the bubble near the original bubble po- 
sition. If, for example, the large circular radio 
emission regions are located at a distance of 50 
kpc from the center (horizontal dotted line) then 
a bubble could rise from 5 kpc to 50 kpc with the 
amount of "entrained" gas corresponding to ap- 
proximately 35% of its initial volume. 



two black curves show, respectively, the boundary 
of the bubble at the beginning of the simulation 
and the region in the cluster where the cluster gas 
has the same entropy as the gas initially inside 
the bubble. It is clear that at the end of the sim- 
ulation, tracer particles are concentrated along a 
surface where cluster gas has a similar entropy to 
that of the gas in the original bubble. Pancake- 
like structures are already visible after ~ 10^ years 
from the start of the simulations. Oscillations near 
the neutral buoyancy point were also observed. 

Strongly underdense bubbles may rise to much 
larger distances from the cluster center before 
reaching a neutral buoyancy point (see Fig. 10). 
We did not simulate the evolution of strongly un- 
derdense bubbles to such late stages, but qualita- 
tively the picture should be the same. 

The largest structures visible in Figure 1 are the 
two almost circular low surface brightness emis- 
sion regions to the northeast and southwest of the 
centre. Owen et al (2000) suggest that these are 
spherical bubbles of radio plasma seen in projec- 
tion. In this case, they may receive substantial en- 
ergy input from the inner lobes which contain the 
active jet. Alternatively one can identify the circu- 
lar emission regions seen in the radio images with 
earlier bubbles that have reached their isodensity 
distance and are now in the process of spreading 
and forming pancake-like structures. We stress 
again that this requires that either rclativistic and 
thermal particles are mixed on microscopic scales 
or the lumps of rclativistic plasma and the ambi- 
ent medium are bound together. 

Owen et al. (2000) have noted the sharp bound- 
ary of the M87 halo which is traced by all observa- 
tions at different frequencies and they argue that 
this excludes a dominant diffusive transport pro- 
cess for the rclativistic plasma. In the proposed 
model the radius and sharpness of the boundary 
is naturally explained by the height to which buoy- 
ant bubbles can rise. Also the complicated struc- 
ture of the outer halo may in this model reflect 
several layers of partly fragmented bubbles. 

The time required for pancake formation at 
large distances from the cluster center is at least 
several times longer than that required for a 
strongly buoyant bubble to travel such distances. 
Therefore, it is extremely difficult for rclativistic 
particles to survive (see 3.3). Either reaccelera- 
tion or exposure of relativistic particles to stronger 
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magnetic fields seems to be inevitable to explain 
the radio emission from the outer radio halo if it is 
indeed the result of lateral spreading of a bubble 
which has reached the neutral buoyancy point. 

If, on the contrary, lumps of radio emitting 
plasma are not bound to the captured lumps of 
thermal gas, then radio bubbles will continue to 
rise along the pressure gradient and will eventu- 
ally leave the central cooling flow region. Over- 
dense lumps of captured gas will separate from 
the bubbles and fall back towards the central re- 
gion. One possible example of overdense material 
having fallen back toward the central galaxy may 
be found in A1795. Cowie et al. (1983) reported 
a remarkable ~ 45" (41 kpc long for the Hubble 
constant of 100 km s^^ Mpc"-"-) filamentary emis- 
sion line (Ha + [NH]) chain originating near the 
central cD galaxy. 

3.6. Total energy budget 

The X-ray luminosity of the whole cooling flow 
region around M87 is about 10^^ erg s~^. On the 
other hand, estimates of the total energy supplied 
by the jet to the ICM give a much larger value 
of 10** erg (e.g. Owen, Eilek and Kassim 
2000). If all this energy can be dissipated into the 
thermal gas in the cooling flow region then the 
cooling flow should instead be a "heating" outflow 
(Owen, Eilek and Kassim 2000), at least at the 
present moment. The situation may be transient 
so that inflow and outflow periods alternate de- 
pending on the strength of the AGN activity (e.g. 
Tucker and David 1997, Binney 1999, Owen, Eilek 
and Kassim 2000). We speculate below on the con- 
sequences of a large input of power in the form of 
relativistic plasma into the cooling flow region. 

Initially the radio lobe, inflated by the jet, ex- 
pands supersonically (e.g. Heinz, Reynolds, Begel- 
man 1998). Later, as the size of the radio lobe 
grows, the expansion becomes subsonic. As we 
discussed above (Section 3.2), a large and strongly 
underdense bubble would rise with approximately 
the Keplerian velocity of ~ 400 km s~^. Assuming 
that this velocity limits the growth of the origi- 
nal bubble due to the jet input power L, one can 
write: p^-j-P()47r';-^?; = L, where Pq is the ambi- 
ent gas pressure, F = 4/3 is the adiabatic index 
of the relativistic gas in the bubble, r is the size 
of the bubble. Substituting L ~ 10** erg s~^, 
Po ~ 2 X 10~^° erg cm~^ and v 400 km s~^ we 



find r ^ 5 kpc, which is approximately the size of 
the inner lobes surrounding the jets. If the sub- 
sonic expansion phase is longer than the duration 
of the initial strongly supersonic expansion, then 
most of the energy deposited by the AGN will go 
into the enthalpy of the expanded lobe and not 
into shock heating the cluster gas. 

The inner lobes break loose (during a time 
f ~ r/?; ~ 10^ years) and form bubbles which then 
subsonically rise through the cluster atmosphere 
(Gull and Northover 1973). At any given moment, 
thermal gas tries to settle in the potential of the 
cluster in order to form a nondecreasing entropy 
profile. Rising bubbles produce convective mo- 
tions in the gas and mix lumps of thermal gas with 
different entropies. Using tracer particles added to 
the thermal gas (see Fig. 9) we estimated that some 
3-4xlO^M0 of cool, thermal gas has been uplifted 
from below 15 kpc to above 35 kpc by the bubble 
after ~ 8 x lO'^ years. Much larger masses of gas 
undergo smaller displacements by that time. For 
example, ~ 1.3x 1O^M0 has been lifted from below 
14 kpc to above 25 kpc. Although these are very 
crude estimates, which are sensitive to the particu- 
lar initial and boundary conditions adopted in our 
simulations, one can expect that if the bubbles are 
launched approximately every 10*" years then they 
are able to transport some 1-10 Mq of cold gas 
per year from the center to the periphery of the 
cooling flow. Interestingly this value is compara- 
ble to the mass deposition rate in M87 estimated 
from the X-ray data. One can further list several 
important consequences of the convection induced 
by the rising bubbles: 

• Due to uplifted cold gas, the X-ray surface 
brightness distribution should be shallower 
than in the absence of bubbles. The entropy 
profiles also flatten. 

• If the bubbles are formed and are rising 
systematically over some preferred direction 
then the cold gas will be continuously up- 
lifted over this direction while it will be flow- 
ing inward in other directions. 

• Thermal instabilities can effectively take 

place within an environment of extensive 
convective motions. As a result, distributed 
mass deposition from cooling gas can natu- 
rally occur. 
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• Ovcrdcnsc lumps of thermal gas, uplifted by 
the rising bubbles, may fall back (when sep- 
arated from the "carrier" buoyant bubble) 
to the central region thereby producing fila- 
mentary structures, stretching radially from 
the galaxy. 

Even if the bubbles are not mixed on micro- 
scopic scales with the thermal plasma and are able 
to eventually leave the cooling flow region, most 
of the energy has to be deposited in the thermal 
gas (if there is a strong contrast in pressure be- 
tween the inner and outer regions of the cooling 
flow). For a bubble moving at a constant ter- 
minal velocity the dissipation of energy is char- 
acterized by the change of the bubble enthalpy 
Y^PV oc {P/Po)~c^, where P is the pressure 
of the cluster thermal gas. The bubble energy 
goes into sound waves, internal waves, turbulent 
motion in the wake, potential energy of the up- 
lifted cold gas etc. Sound waves can carry energy 
away from the central region, while most of the 
other channels would eventually result in heating 
the cooling flow region. If the jet power is indeed 
^ 10^^ erg s"-'^, then even a modest 10% efficiency 
of local dissipation into heat should be enough 
to exceed the radiative cooling of the gas. The 
gas then heats up, but this heating is distributed 
through a large volume. Since no shock waves are 
envisaged in this picture one would not expect to 
find localized, hot, compressed regions. 

If strong energy input from AGN is maintained 
for a very long period (and dissipation into heat 
is efficient), then a core with a flat entropy profile 
will be formed and the size of this core will grow 
with time. We note here that the rate of energy 
extraction from the bubble depends on the pres- 
sure gradient (see above expression for enthalpy). 
Therefore, if the gas is heated to a temperature 
comparable to the depth of the local potential well, 
the dissipation efficiency will decrease in this re- 
gion. This effect may provide a self-regulation for 
heating in the central part of the cooling flow re- 
gion. 

Eventually, strong heating is capable of com- 
pletely eliminating the cooling flow, unless the 
power of an AGN is in turn regulated by the con- 
ditions in the cooling flow region as discussed by 
Tucker and David (1997), Binney (1999), Owen, 
Eilek and Kassim (2000). In M87, heating at a 



rate of ~ 5 x 10'*^ ergs would eliminate the 
cooling flow for ~ 10^ years. This is compara- 
ble to the age of the largest radio structures. On 
the other hand, the complicated radio halo mor- 
phology suggests that periods of high AGN activ- 
ity have alternated with the periods of low AGN 
power. We note here that dissipation of the sub- 
sonic motions into heat may take a longer time 
than the characteristic bubble rise time and there- 
fore the actual gas heating rate may be related to 
the AGN activity averaged over a long period of 
time. 

4. Conclusions 

Buoyant bubbles of cosmic rays, slowly rising 
through the cooling gas, can qualitatively explain 
the complicated radio and X-ray morphology of 
the central 40 kpc region around M87. Torus- 
like features, seen in the radio, may be similar 
to the "mushrooms" which appear in Rayleigh- 
Taylor unstable conflgurations: as the fluid rises 
through the ambient medium, Kelvin-Helmholtz 
instabilities create torus-like "mushroom" heads. 
The excess X-ray emission trailing some promi- 
nent radio features could be due to cold gas cap- 
tured by the rising bubbles and uplifted to large 
distances from the central source. 

Rising bubbles produce strong convection in 
the cooling flow region. Convection flattens X- 
ray surface brightness and entropy profiles in the 
cooling flow. Convection may promote distributed 
mass deposition in the cooling flow and the for- 
mation of the filaments. If more than 10% of the 
energy supplied by the jet is dissipated into heat 
in the cooling fiow region, then large cores with 
flat entropy profiles may be formed. 

New Chandra and XMM observations of cool- 
ing flow clusters (e.g. McNamara et al. 2000, 
Fabian ct al. 2000a) demonstrate that the in- 
teraction of the radio sources and thermal gas 
is widespread in these objects. Furthermore, the 
mass of cooling gas at temperatures less than ^1 
keV, derived from the spectra, seems to be sig- 
nificantly smaller than expected from the sim- 
plest cooling flow models (Bohringer et al. 2001, 
Tamura et al. 2001, Peterson et al. 2001). Pos- 
sible explanations for this discrepancy are sum- 
marized by Fabian et al. (2000b; sec also refer- 
ences therein). As we discussed above, heating 
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Fig. 11. — Simple simulations demonstrating formation of a pancake. The top panel shows the distribution 
of the quantity T jr?^"^ which characterizes the entropy of the gas. Bottom panel shows the distribution of 
tracer particles, originally located within the bubble. The left column corresponds to the state shortly after 
the beginning of simulations and the right column corresponds to the state after ~ 2.6 x 10* years. The 
initial entropy of the bubble (green circle in the top-left plot) was such that the cluster gas has a similar 
entropy at a radius of 25 kpc from the center. The entropy contrast between the bubble and surrounding 
cluster gas is rather modest and the bubble travels only a small distance from the center. The contours of 
the initial bubble and a circle with a 25 kpc radius are shown by black lines. The outer edge of the simulated 
volume (outer black line in the bottom panel) corresponds to a distance of 50 kpc from the cluster center. 
It is clear that at the end of the simulation, tracer particles are concentrated along the surface where cluster 
gas has a similar entropy to that of the gas in the original bubble. 



by the buoyant bubbles may play an important 
role in the thermal balance in cooling flows. One 
of the predictions of the buoyant bubble picture, 
namely a correlated morphology of the colder up- 
lifted gas and radio emission, seems to be in a 
broad agreement with recent XMM results (Bel- 
sole et al. 2001). 

Our model is of course highly oversimplified 
and is not intended to closely reproduce the M87 
environment. Clearly more detailed simulations 
are needed, perhaps in 3D, which would also fol- 
low the formation of the initial bubble (possibly 
with continuous energy input) and include radia- 
tive cooling. However, despite the simple treat- 
ment of bubbles in cluster and galaxy atmospheres 
presented above, we believe that this approach in- 
dicates areas of particular interest. 

We are grateful to the editor Steven N. Shore 
and the referee for important comments and sug- 
gestions. We thank Alcxei Kritsuk, Nail In- 



ogamov, Christine Jones, Nail SibgatuUin and 
Henk Spruit for useful discussions. This re- 
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the High Energy Astrophysics Science Archive Re- 
search Center Online Service, provided by the 
NASA/Goddard Space Fhght Center. 
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